%**************************************************************************
%
%           "Endogenous Liquidity and Capital Reallocation"
%      Cui, W., Wright, R., & Zhu, Y. (2024), Journal of Political Economy
%                       Last Modified: Feb 2024.  
%
%**************************************************************************

%%%% Computing steady state equilibrium
%%%% I.I.D. Model case with credit and money

clc; clear; close all; format short


options = optimset('Display','off'); 
options.OptimalityTolerance = 1e-8; 
options.StepTolerance       = 1e-8;
options.MaxIter             = 50000; 
options.MaxFunEvals         = 50000; 
options.FunctionTolerance   = 1e-8;


%% 1. ---------------------   Initilization  ------------------------------
par.credit_parameter= 'CHI_PI'; % calibrate a particular credit parameter and the rest are zeros
par.Rel_Tol         = 1e-5; % for numerical integration
par.Abs_Tol         = 1e-6; % for numerical integration
par.int_method      = 'auto'; % for numerical integration

%%%% fixed parameters across models (some are initial guess, to be calibrated)
par.BETA            = 0.9619;   % discount factor
par.SIGMA           = 1.0000;   % CRRA parameter
par.ETA             = 0.6110;   % labor share (1 - ETA calibrated to investment/output ratio)
par.DELTA           = 0.1000;   % depreciation rate; (implied capital share to target targeting I/Y ratio 0.20; mannually adjust)
par.A               = 1;        % normalized aggregate productivity
par.tau_k           = 0.2500;   % capital tax rate
par.tau_h           = 0.2200 ;  % labor tax
par.iota            = 0.0672;   % interest rate (steady state)
par.ALPHA           = 1;        % matching probability
par.eps_            = 0.0001;   % the deviation for calculating AQC spending derivatives cannot be too large (e.g., <0.01)
par.iota_per_change = par.eps_ / (1 + par.iota);

par.log_eps_mu      = 0;        % log-normal capital productivity mean
par.log_eps_sig     = 0.5 / (1 - par.ETA); % log-normal capital productivity std.dev
par.MU              = (1 + par.iota) * par.BETA - 1;  % the corresponding steady-state money supply/inflation rate
par.S_PARAM         = par.log_eps_sig;
par.M_PARAM         = par.log_eps_mu;
par.EPS_L_BOUND     = 0.0001;
par.EPS_H_BOUND     = 10000;
par.a               = 0;        % IID case; productivity of the persistent part
par.THETA           = 0.50;     % bargaining power

%%%% targets, and INITIAL values for parameters to hit targets
% Note: Y always stands for output without financial services!
par.H_target        = 1 / 3;
par.R_ratio_target  = 0.322;
par.P_share_target  = 0.24;
par.IM_share        = 0.011;  % intermediation share of the output (that does not contain intermediation)
par.Cash_target     = 0.042 * (1 + par.IM_share);     % cash to (CM) output ratio; cash to total output (excluding imports/exports) is 4.2% in the data
par.G_Y             = 0.160 * (1 + par.IM_share);     % government spending to (CM) output ratio; G over total output = 0.16 in the data
par.aqc_elasticity_target = -0.64;

% the following parameters are "initial values"; they will be calculated
par.XI              = 2.3000;   % disutility parameter for working; initial guess 
par.CHI_0           = 0;        % non-secured credit parameter
par.CHI_k           = 0;        % existing capital collaterability
par.CHI_PI          = 0.1;      % Holmostrom-Tirole credit parameter
par.CHI_q           = 0.90;     % collaterability: WLOG as we have money as substitutes
par.c_mu            = -1.5;     % cost of entry mean, \mu_{\gamma} in the paper notation
par.c_sig           = 2;        % cost of entry std.dev., \mu_{\sigma} in the paper notation


%% 2. --- Define distribution and functions with idiosyncratic shocks -----

fun                 = define_function_fun(par);
par.q_average_prod  = integral(@(x) x .* fun.f(x, par.a), fun.l_bound(par.a), fun.h_bound(par.a)); % the average productivity

upper_bound_mean    = [integral(@(x) fun.f(x, par.a), fun.l_bound(par.a), fun.h_bound(par.a)),...
integral(@(x) x.* fun.f(x, par.a), fun.l_bound(par.a), fun.h_bound(par.a))];

fprintf('Pre-analysis Check: CDF at upper bound and the average productivity %.3f %.3f', upper_bound_mean)
disp(' '); disp(' ')
fprintf('Note: Theoretically, they should be 1 and %.3f', exp(par.log_eps_mu + 0.5 * par.log_eps_sig^2));
disp(' '); disp(' ')


%% 3. ----------    Run CWZ model, steady state      ----------------------

disp('******************>Main Model Calibration <************************')
disp(' ')

%%%% produce calibrated parameters stored also in par
do_calibration_main_model;

%if THETA == 0.5
%save param par

eqm_calibrated = eqm;
par_calibrated = par;

save eqm_calibrated eqm_calibrated
save param_calibrated par_calibrated

%end


